package PlotPourri;

import PlotPourri.Database.DatabaseSMCi;
import PlotPourri.Database.StarDatabase;
import com.panayotis.gnuplot.terminal.ImageTerminal;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileNotFoundException;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math3.stat.regression.SimpleRegression;

public class Multiphase {
    
    static int phases = 100;

    public static void main(String[] args) {
        //StarDatabase starDatabase = new StarDatabase(new File("/home/earl/Documents/Research/SMC/TypeI/SMCOGLEIII.dat"));
        //StarDatabase starDatabase = new LMCiDatabase(new File("E:\\Research\\Astrophysics\\LMC\\TypeI\\OG3.dat"));
        StarDatabase starDatabase = new DatabaseSMCi(new File("E:\\Research\\Astrophysics\\SMC\\TypeI\\SMCOGLEIII.dat"));
        
        int order = 8;

        List<List> plots = new ArrayList<List>();
        for (int i = 0; i <= phases - 1; i++) plots.add(new ArrayList<double[]>());
        for (int name : starDatabase.keySet()) {
            Star star = starDatabase.get(name);
            double log10P = Math.log10(star.getPeriod());
            if (log10P < /*0.4*/1) continue;
//            double I = star.Ao(order, BandName.I);
//            double V = star.Ao(order, BandName.V);
//            System.out.println(star.getPeriod() + " " + (I - 1.55 * (V - I)));
            double[][] I = star.getBand(BandName.I).get(0).getShiftedData(order, phases);
            double[][] V = star.getBand(BandName.V).get(0).getShiftedData(order, phases);
            //double[][] data = V;
            //for (int i = 0; i < I.length - 1; i++) plots.get(i).add(new double[] {log10P, V[i][1] - I[i][1]}); // Period-Color
            //for (int i = 0; i < data.length - 1; i++)
            //    if (data[i][0] >= 0 && data[i][0] <= 1 && data[i][1] >= 10 && data[i][1] <= 20)
            //        plots.get(i).add(new double[] {log10P, data[i][1]}); // Period-Luminosity
            for (int i = 0; i < I.length - 1; i++)
                if (I[i][0] >= 0 && I[i][0] <= 1 && I[i][1] >= 10 && I[i][1] <= 20 && V[i][0] >= 0 && V[i][0] <= 1 && V[i][1] >= 10 && V[i][1] <= 20)
                    plots.get(i).add(new double[] {log10P, I[i][1] - 1.55 * (V[i][1] - I[i][1])}); // Wessenheit
        }
        
        for (int i = 0; i < phases; i++) for (int j = i + 1; j < phases; j++) System.out.println(i + " " + j + " " + PearsonCC(plots.get(i), plots.get(j)));

//        List allList = new ArrayList<Double>();
//        List leftList = new ArrayList<Double>();
//        List rightList = new ArrayList<Double>();
//        List allError = new ArrayList<Double>();
//        List leftError = new ArrayList<Double>();
//        List rightError = new ArrayList<Double>();
//        System.out.printf("# %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \n", "i", "<10 Slope", "Error", "<10 Intercept", "Error", ">10 Slope", "Error", ">10 Intercept", "Error");
//        for (int i = 0; i < plots.size(); i++) {
//            List<double[]> plot = plots.get(i);
//            double[][] data = plot.toArray(new double[plot.size()][2]);
//            SimpleRegression all = new SimpleRegression();
//            SimpleRegression left = new SimpleRegression();
//            SimpleRegression right = new SimpleRegression();
//            for (double[] point : data) {
//                if (point[0] > 0.4 && point[0] < 1) left.addData(point[0]-0.4, point[1]);
//                else if (point[0] >= 1) right.addData(point[0]-1, point[1]);
//                all.addData(point[0], point[1]);
//            }
//            allList.add(all.getSlope());
//            leftList.add(left.getSlope());
//            rightList.add(right.getSlope());
//            allError.add(all.getSlopeStdErr());
//            leftError.add(left.getSlopeStdErr());
//            rightError.add(right.getSlopeStdErr());
//            System.out.printf("%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \n", i / (double) 100, left.getSlope(), left.getSlopeStdErr(), left.getIntercept(), left.getInterceptStdErr(), right.getSlope(), right.getSlopeStdErr(), right.getIntercept(), right.getInterceptStdErr());
//        }
//        System.out.println("all data: " + average(allList) + " +- " + average(allError));
//        System.out.println("<~6.30957days: " + average(leftList) + " +- " + average(leftError));
//        System.out.println(">~6.30957days: " + average(rightList) + " +- " + average(rightError));

        //JavaPlot p = new JavaPlot();
        //printData(plots);
/*        for (int i = 0; i < plots.size(); i++) {
            JavaPlot p = new JavaPlot();
            List<double[]> plot = plots.get(i);
            double[][] data = plot.toArray(new double[plot.size()][2]);

            DataSetPlot dsp = new DataSetPlot(data);
            dsp.setTitle("Data");
            p.addPlot(dsp);

            SimpleRegression left = new SimpleRegression();
            SimpleRegression right = new SimpleRegression();
            SimpleRegression all = new SimpleRegression();

            for (double[] point : data) {
                if (point[0] < 1) left.addData(point[0], point[1]);
                else right.addData(point[0], point[1]);
                all.addData(point[0], point[1]);
            }

            FunctionPlot line = new FunctionPlot(all.getSlope() + "*x+" + all.getIntercept());
            line.setTitle("Linear Regression");
            line.set("lw", "3");
            p.addPlot(line);

            FunctionPlot leftLine = new FunctionPlot("(x < 1 ? " + left.getSlope() + "*x+" + left.getIntercept() + " : 1/0)");
            FunctionPlot rightLine = new FunctionPlot("(x >= 1 ? " + right.getSlope() + "*x+" + right.getIntercept() + " : 1/0)");
            leftLine.setTitle("< 10 Days");
            rightLine.setTitle("> 10 Days");
            leftLine.set("lw", "3");
            rightLine.set("lw", "3");
            p.addPlot(leftLine);
            p.addPlot(rightLine);

            //FunctionPlot piecewise = new FunctionPlot("(x < 1 ? " + left.getSlope() + "*x+" + left.getIntercept() + " : " + right.getSlope() + "*x+" + right.getIntercept() + ")");
            //p.addPlot(piecewise);

            String phase = i + "";
            while (phase.length() < 3) phase = "0" + phase;
            File file = new File("/home/earl/Documents/Research/SMC/TypeI/W-data/" + phase + ".png");
            ImageTerminal png = new ImageTerminal();
            createFile(file, png);

            // Period-Luminosity
//            p.setTitle("Period-Luminosity Relation for V Band at Phase " + phase);
//            p.set("ylabel", "\"mag(" + phase + ")\"");
//            p.set("xlabel", "\"log P\"");
//            p.set("yrange", "[18:8]");
//            p.set("xrange", "[0.4:2]");

            // Period-Color
//            p.setTitle("Period-Color Relation at Phase " + name);
//            p.set("ylabel", "\"(V(" + name + ") - I(" + name + "))\"");
//            p.set("xlabel", "\"log P\"");
//            p.set("yrange", "[0:2]");
//            p.set("xrange", "[0.4:2]");

            p.setTitle("Period-Wesenheit Relation in the SMC at Phase " + phase);
            p.set("ylabel", "\"I - 1.55 * (V - I)\"");
            p.set("xlabel", "\"log P\"");
            p.set("yrange", "[18:8]");
            p.set("xrange", "[0.4:2]");
            
            p.set("key", "right bottom");
            p.setTerminal(png);
            p.plot();
            try {
                ImageIO.write(png.getImage(), "png", file);
            } catch (IOException ex) {
                System.err.print(ex);
            }
            p = null;
        }*/
    }
    
    public static double PearsonCC(List<double[]> X, List<double[]> Y) {
        double SumXi = 0, SumXiSq = 0, SumYi = 0, SumYiSq = 0, SumXiYi = 0;
        for (int i = 0; i < phases; i++) {
            SumXi += X.get(i)[1];
            SumXiSq += Math.pow(X.get(i)[1], 2);
            SumYi += Y.get(i)[1];
            SumYiSq += Math.pow(Y.get(i)[1], 2);
            SumXiYi += X.get(i)[1] * Y.get(i)[1];
        }
        System.out.println("SumXi = " + SumXi);
        System.out.println("SumXiSq = " + SumXiSq);
        System.out.println("SumYi = " + SumYi);
        System.out.println("SumYiSq = " + SumYiSq);
        System.out.println("SumXiYi = " + SumXiYi);
        return (phases * SumXiYi - SumXi * SumYi) / (Math.sqrt(phases * SumXiSq - Math.pow(SumXi, 2)) * Math.sqrt(phases * SumYiSq - Math.pow(SumYi, 2)));
    }

    static double average(List<Double> list) {
        double sum = 0;
        for (Double d : list) sum += d;
        return (list.size() > 0) ? sum/(double) list.size() : 0;
    }

    static void printData(List<List> plots) {
        System.out.printf("# %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \t %s \n", "i", "<10 Slope", "Error", "<10 Intercept", "Error", ">10 Slope", "Error", ">10 Intercept", "Error");
        for (int i = 0; i < plots.size(); i++) {
            List<double[]> plot = plots.get(i);
            double[][] data = plot.toArray(new double[plot.size()][2]);
            SimpleRegression left = new SimpleRegression();
            SimpleRegression right = new SimpleRegression();
            for (double[] point : data)
                if (point[0] > 0.4 && point[0] < 1) left.addData(point[0], point[1]);
                else if (point[0] > 1) right.addData(point[0], point[1]);
            System.out.printf("%f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \t %f \n", i / (double) 100, left.getSlope(), left.getSlopeStdErr(), left.getIntercept(), left.getInterceptStdErr(), right.getSlope(), right.getSlopeStdErr(), right.getIntercept(), right.getInterceptStdErr());
        }
    }

    static void createFile(File file, ImageTerminal terminal) {
        try {
            file.createNewFile();
            terminal.processOutput(new FileInputStream(file));
        } catch (FileNotFoundException ex) {
            System.err.print(ex);
        } catch (IOException ex) {
            System.err.print(ex);
        }
    }
}

